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Abstract 

We give a general method to calculate photonic band structure in the form of wave 
number k as a function of frequency u, which is required whenever we want to calculate 
c/3 ' signal intensity related with photonic band structure. This method is based on the fact 
that the elements of the coefficient matrix for the plane wave expansion of the Maxwell 
equations contain wave number up to the second order, which allows us to rewrite the 
original eigenvalue equation for uj"^ into that for wave number. This method is much 
I better, especially for complex wave numbers, than the transfer matrix method of Pendry, 
which gives the eigenvalues in the form of exp[ikd] . As a simplest example, we show a 
^ ■ comparison of uj{k) and k{uj) for a model of intersecting square rods. 

O 

; 1 Introduction 

^ ■ Usual band structure is calculated as eigen frequency for a given wave vector of a photonic 
crystal. This is common to the electronic band structure of crystals, too. The information 
about {hujj{k)} serves as a fundamental quantity to discuss the behavior of the crystal. 
Q I However, when we want to calculate any signals from this crystal, we have to work on 
a sample with boundary, and consider the relationship between the solutions in- and 
outside the sample. For both photonic and electronic bands, we need to connect the 
solutions across the boundary, and in doing this we have to consider both propagating 
[ and evanescent wave components. For this purpose, we have to find out a method to solve 
Q I the eigenvalue equation, not in the form of Wj{k), but in the form of kj{uj) for real u. 
J-> ■ The eigenvalue equation can be expressed in general as a linear equations for the 

Fourier components {X} of wave functions or EM field 

SX = (1) 

for both electronic and photonic band problems. This is the Schrodinger equation for 
electronic problems, and Maxwell equations for photonic crystals. The coefficient matrix 
is a function of {uj,k) as well as Fourier components of "potential" term, which depend 
on the reciprocal lattice vectors. This matrix has the form 

S = S' -XI (2) 



1 



where X = hcu (energy) for Schrodinger equation, and A = cj^ for Maxwell equations. 
Thus the eigenvalues of the matrix S' give the band structure {ujj{k)} in a straightforward 
manner. 

On the other hand, it is not so clear to find a way to get kj{uj) from the same equation. 
However, if we note the following simple argument, which is known for a numerical solution 
of eigenvalue equation, we are lead to an appropriate way to get kjiuj). 

Let us consider a set of equations 

A X = . (3) 

where the matrix A is assumed to depend on a parameter A. We now ask the condition 
for A to obtain a nontrivial solution X. If the A-dependence of the matrix A is 

A = B - A C , (4) 

the condition for the nontrivial solution, i.e., the eigenvalue equation for A, is written as 

det|C"^ B- Al| = . (5) 

In a similar way, for a matrix A given in the following form, 

A = B + A C - A^ D , (6) 

we rewrite eq.(3) by using new variables Y — XX as 

D-^BX + D-^cy - Ay = (7) 

Then, the condition to have nontrivial solution for {X, Y}, which works as an eigenvalue 
equation for A, can be written as det|U — Al| = , where 



U 



, 1 
D-^B , D-^C 



(8) 



This way of rewriting is possible for any power of A, leading to the eigenvalue problem 
for A. Being well-known as a technique to get appropriate form of eigenvalue equation, it 
has been used for the calculation of multibranch polariton dispersion curves in the form 

0ik{Lj) [1]. 

The cases of electronic and photonic band structure calculation in plane wave expan- 
sion have A (= A;) up to the second order. This allows us to write down the general 
matrix equation for the arbitrary configuration as shown in detail below, where we will 
treat the case of photonic bands explicitly. Its application to the case of electronic bands 
is straightforward, as mentioned above. 
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2 Formulation 



If we eliminate the electric field in the Maxwell equations for a photonic crystal ("H- 
method" [2]), and make a Fourier series expansion with respect to the reciprocal lattice 
vectors G of the crystal, we get a set of equations 

v G' 

where 77) are the Cartesian coordinates, and 

Sg,G' = ^Vg-G' + (Ve)G-G' (fc + G) X (fc + G') x . (10) 

In these expressions, e and ^ are the dielectric constant and magnetic permeability with 
periodic structure of the photonic crystal in consideration, and qi the magnetic field 
of light in Fourier representation. For a perfect crystal, fc is a good quantum number, 
and det|S| = gives the nontrivial solutions, i.e., the dispersion relation {ijjj{k)} of the 
eigenmodes. 

It should be noted that the coefficient matrix S contains k up to the second order. 
Let us introduce a surface of this photonic crystal, which is always necessary to define 
an observation process. For simplicity we assume that the outside is vacuum. We take 
the z axis to the surface normal direction (toward inside) , and the incident light is within 
the xz plane. Then, the incident wave vector has components (/cy, 0, k). If we denote the 
reciprocal lattice vectors of the surface periodic 2D lattice as {g}, there are reflected and 
transmitted light beams with the lateral wave vector components of {fc|| + g}. The z- 
components of the wave vectors of reflected and transmitted lights should be determined 
so as they satisfy the dispersion relations in vacuum and bulk crystal, respectively. An 
important feature about the 2;-components is that they may well be complex numbers, 
namely, they can be evanescent waves. Because of the presence of the surface, which 
breaks the translational symmetry, the evanescent solutions are allowed as long as their 
amplitudes decay as they go away from the surface. 

Noting the identity 

(fe+G)x(fc+G')xi?fc_(., = {k+G') {{k+GyHj^c^')-Hk,G' {{k+Gy{k+G')] (11) 
and 

fc + G = {k\\ + G^)e^ + GyQy + (/« + GJe^, , (12) 
we obtain the components of S for a given block of (G, G') as 

S = '^Vg-G'I + {-(^ + G)-{k + G')l + Si}(l/6)^_Gf' , (13) 

where 



Si 



m + G',)m + G.) , (A-ii + G',)Gy , (A;|| + G',){k + 

G'y{k\\+G^), G'yGy, G'y{K + G, 

{K + G'^)m + G,) , («: + G;)G, , («: + G'J(«: + G, 



(14) 
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Equation (9) should be supplemented by the condition V • .ff = 0, which requires the 
transversal nature of the magnetic field at any point. In terms of the Fourier components, 
this condition leads, for all G's, to 

ih + G'^)H^^]^, + G'X^]^' + {^ + G'.)Hi]G' = 0. (15) 

These equations are used to ehminate {H^^^qi} from eq.(9), which leads to the final form 
of the equations to be solved as 

If we take N reciprocal lattice vectors, we need to determine 2N components {Hj^^^,, H^^^,} 
Choosing y and z for ^, wc explicitly write this matrix equation in the form of eq.(6), 
which can be transformed into the eigenvalue problem for k as discussed above. The size 
of the matrix equation for k is 4iV x 4N, which gives 2N solutions for each of the positive 
and negative directions of k, including propagating and evanescent modes. The factor 2 
corresponds to the polarization direction. 



Fig.l shows an example of comparison of the band calculations according to this 
method and H-method for a photonic bands of intersecting square rods for propagating 
modes. The number of the reciprocal lattice vectors used for this calculation is N — 
9x9x9. As expected, the calculated bands as uj{k) and k{u) agrees very well. Fig. 2 
shows the bands of complex wave vectors. The branches with large imaginary part are 
obtained as smooth curves, indicating the stability of numerical calculation. However, 
the fact that the lm[K] can exceed 27r x 6 for the present choice of N indicates that in 
terms of the transfer matrix method we have to deal with the eigenvalues ranging from 
10^^ to 10~^^. For higher frequency regions, we need shorter wavelength components, i.e., 
more reciprocal lattice vectors, which gives steeper evanescent wave components. Thus 
the difficulty to solve eigenvalue problem of k increases. In this way, the practical limit 
is reached rather quickly if we use the transfer matrix method. The reason why we need 
more reciprocal lattice vectors for higher frequency is that the boundary conditions for a 
higher frequency light should be given at finer spatial mesh points, because it has shorter 
wavelength. 



3 Discussion 

3.1 Comparison with transfer matrix method 

It is well known that the band calculation to get {kj{uj)} can be done via so-called 
transfer matrix method [3]. In this method, one discretizes the Maxwell equations in 
an appropriate mesh, which provide the relationship among the EM field components 
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ascribed to these mesh points. From this set of difference equations, one can derive a 
transfer matrix T which relates the EM fields at various positions on one lattice plane 
with those on the neighboring lattice plane. Since the EM fields on the corresponding 
positions on the neighboring lattice planes (with a spacing d) are related according to 
Bloch theorem, the eigenvalues of T should be exp[ife • d], i.e., 

det|T-exp[ife-d]| =0 . (17) 

This logic is certainly correct, but it is not quite convenient to carry it out as a 
numerical calculation. Since we need all the possible solutions of the band states including 
evanescent waves, the eigenvalues exp[ife-<i] can be extremely large or small, which causes a 
serious trouble to the numerical calculation. This means that the accuracy of calculating k 
can become rather poor. This has been encountered in our previous work to determine the 
branching ratio of an incident light into various waves in the bulk [5] . For lower frequency 
region, the result is not too bad. However, as we examine higher frequency region, we 
start to see the frequencies where the reflectivity exceeds unity, i.e., a typical violation of 
physics principle. We beheve that this is caused by the inaccuracy of calculating fc's for 
a given frequency. 

In contrast, the formulation in the previous section gives an eigenvalue equation for 
k{u!), not for exp[ifi;ci]. Therefore, the range of eigenvalues is not tremendously large, as 
shown in Fig. 2, and we can hope a good accuracy in their numerical values. 

3.2 Branching ratio in the reflected and transmitted light beams 

For a given incident light, there arise various beams of reflected and transmitted lights with 
propagating and evanescent character due to (i) the multibranch structure of photonic 
bands and (ii) the Bragg scattering from the surface periodicity. In order to fix the 
amplitude ratio of these different beams, it is enough to consider the Maxwell boundary 
conditions (MBC's) at the surface. In spite of the multibranch structure of photonic 
bands, there is no need of "additional boundary condition (ABC)" [4], because we are 
dealing with macroscopic local response without any spatial dispersion effect. 

Since the surface is not uniform, we need to require MBC's at each point. Alternatively, 
we may require the MBC's for each Fourier components of electric and magnetic fields 
with respect to the surface reciprocal lattice vectors {g}. The number of gi's to be used 
should be determined by the desired accuracy of numerical calculation. For each g, we 
take the corresponding beams on both sides of the surface, and make them fulfil the 
MBC's. In this way we can determine the amplitudes of all the beams uniquely. 

It should be noted that any beam in the photonic crystal has the form of Bloch 
wave, i.e., a product of fc-dependent phase factor and lattice periodic part. Therefore, 
its amplitude at the surface changes according to the position of a surface cross section 
for any given Miller indices. This means that the beam amplitudes of transmitted and 
reflected lights depend on the surface cross section for any given Miller indices. This was 
explicitly demonstrated by us for the [100] surface of the intersecting square rods, which 
shows remarkable difference among different cross sections of the surfaces [5]. This fact 
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shows a possibility of designing the surface to control the flow of light energy in desired 
beams. For this calculation also, the new method of calculating k{tu) should be quite 
useful, because the boundary conditions depend on the accuracy of evanescent waves. 

3.3 Other applications 

In more realistic situations, we have to consider the finite size of photonic crystals, or even 
systems without periodicity. In such cases, the method of FDTD (Finite Difference in 
Time Domain) method [6] is well known and used quite popularly. The present method 
may be used in such cases, if we put the matter in a periodic system with super-cell 
construction. In this way, we can hopefully extend the present scheme to a wide class of 
matter systems to calculate their optical responses in frequency domain in a straightfor- 
ward way. This will provide us with a complimentary information from that of FDTD 
method. 

For the discussion of nonlinear processes in photonic crystals, we also need to know 
how much of an incident light is fed by a particular branch of photonic bands. Thus, we 
need to prepare this input information by the argument of boundary conditions, and this 
is important both for fundamental and applicational purposes. 
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Figure 1: Band structure of an intersecting square-rod model calculated as uj{k) and k{uj) 
for real k = (0.1, 0, K,)27i/d. The number of the reciprocal lattice points is 9 x 9 x 9, and 
other parameters are the same as in ref.[5]. 
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Figure 2: Imaginary parts of the evanescent solutions of fe(a;) for the same frequency 
range as in Fig.l. The number of branches in this region is 2 x 9 x 9, where the factor 2 
counts the polarization degree of freedom, including real (propagating) modes. 
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